#include "Sun.h"

CSun::CSun()
{
	char timebuf[128];
	_strtime_s(timebuf);
	int hour = (timebuf[0] - '0') * 10 + (timebuf[1] - '0');
	int minute = (timebuf[3] - '0') * 10 + (timebuf[4] - '0');
	_strdate_s(timebuf);
	int month = (timebuf[0] - '0') * 10 + (timebuf[1] - '0');
	int day = (timebuf[3] - '0') * 10 + (timebuf[4] - '0');
	int year = (timebuf[6] - '0') * 10 + (timebuf[7] - '0') + 2000;
	m_date.SetDate(8, year, month, day, hour, minute);
	m_geoInfo.m_latitude = 40;
	m_geoInfo.m_longitude = 116;
	CalcSolarPosition();
}

CSun::CSun(CDate date, CGeoInfo geoInfo)
{
	m_date = date;
	m_geoInfo = geoInfo;
	CalcSolarPosition();
}

void CSun::SetDate(CDate date)
{
	m_date = date;
	CalcSolarPosition();
}

void CSun::SetGeoInfo(CGeoInfo geoInfo)
{
	m_geoInfo = geoInfo;
	CalcSolarPosition();
}

CSpherePosition CSun::GetSolarPosition()
{
	return m_solarPosition;
}

CPoint3D CSun::GetGloSolarPosition()
{
	return m_gloPosition;
}

void CSun::CalcSolarPosition()
{
	// calculate julian day
	double julianDay;
	julianDay = Date2Num(m_date) + 2415018.5 + (m_date.m_hour + m_date.m_minute / 60.0) / 24 - m_date.m_timeZone / 24.0;

	// calculate julian century
	double julianCentury;
	julianCentury = (julianDay - 2451545) / 36525;

	// calculate geom mean long sun angle
	double geomMeanLongSunAngle;
	geomMeanLongSunAngle = fmod(280.46646 + julianCentury * (36000.76983 + julianCentury * 0.0003032), 360);

	// calculate geom mean anom sun angle
	double geomMeanAnomSunAngle;
	geomMeanAnomSunAngle = 357.52911 + julianCentury * (35999.05029 - 0.0001537 * julianCentury);

	// calculate sun eq of ctr
	double sunEqOfCtr;
	sunEqOfCtr = sin(DEG_TO_RAD(geomMeanAnomSunAngle)) * (1.914602 - julianCentury * (0.004817 + 0.000014 * julianCentury)) 
		+ sin(2 * DEG_TO_RAD(geomMeanAnomSunAngle)) * (0.019993 - 0.000101 * julianCentury) + sin(3 * DEG_TO_RAD(geomMeanAnomSunAngle)) * 0.000289;

	// calculate sun true long angle
	double sunTrueLongAngle;
	sunTrueLongAngle = geomMeanLongSunAngle + sunEqOfCtr;

	// calculate mean obliq ecliptic angle
	double meanObliqEclipticAngle;
	meanObliqEclipticAngle = 23 + (26 + ((21.448 - julianCentury * 
		(46.815 + julianCentury * (0.00059 - julianCentury * 0.001813)))) / 60) / 60;

	// calculate sun app long angle
	double sunAppLongAngle;
	sunAppLongAngle = sunTrueLongAngle - 0.00569 - 0.00478 * sin(DEG_TO_RAD(125.04 - 1934.136 * julianCentury));

	// calculate obliq corr angle
	double obliqCorrAngle;
	obliqCorrAngle = meanObliqEclipticAngle + 0.00256 * cos(DEG_TO_RAD(125.04 - 1934.136 * julianCentury));

	// calculate sun decline angle
	double sunDeclineAngle;
	sunDeclineAngle = RAD_TO_DEG(asin(sin(DEG_TO_RAD(sunAppLongAngle)) * sin(DEG_TO_RAD(obliqCorrAngle))));

	// calculate eccent earth orbit
	double eccentEarthOrbit;
	eccentEarthOrbit = 0.016708634 - julianCentury * (0.000042037 + 0.0000001267 * julianCentury);

	// calculate var y
	double y = tan(DEG_TO_RAD(obliqCorrAngle) / 2) * tan(DEG_TO_RAD(obliqCorrAngle) / 2);

	// calculate eq of time
	double eqOfTimeMinute;
	eqOfTimeMinute = 4 * RAD_TO_DEG(y * sin(2 * DEG_TO_RAD(geomMeanLongSunAngle)) - 2 * eccentEarthOrbit * sin(DEG_TO_RAD(geomMeanAnomSunAngle)) 
		+ 4 * eccentEarthOrbit * y * sin(DEG_TO_RAD(geomMeanAnomSunAngle)) * cos(2 * DEG_TO_RAD(geomMeanLongSunAngle)) - 0.5 * y * y * sin(4 * DEG_TO_RAD(geomMeanLongSunAngle)) 
		- 1.25 * eccentEarthOrbit * eccentEarthOrbit * sin(2 * DEG_TO_RAD(geomMeanAnomSunAngle)));

	// calculate true solar time
	double trueSolarMinute;
	trueSolarMinute = fmod((m_date.m_hour * 60 + m_date.m_minute)  + eqOfTimeMinute + 4 * m_geoInfo.m_longitude - 60 * m_date.m_timeZone, 1440);

	if (trueSolarMinute < 0)
	{
		trueSolarMinute += 1440;
	}

	// calculate hour angle
	double hourAngle;
	if (trueSolarMinute < 0)
	{
		hourAngle = trueSolarMinute / 4 + 180;
	}
	else
	{
		hourAngle = trueSolarMinute / 4 - 180;
	}

	// calculate zenith angle
	double solarZenithAngle;
	solarZenithAngle = acos(sin(DEG_TO_RAD(m_geoInfo.m_latitude)) * sin(DEG_TO_RAD(sunDeclineAngle)) 
		+ cos(DEG_TO_RAD(m_geoInfo.m_latitude)) * cos(DEG_TO_RAD(sunDeclineAngle)) * cos(DEG_TO_RAD(hourAngle)));

	m_solarPosition.m_elevation = (float)(M_PI / 2 - solarZenithAngle);

	solarZenithAngle = RAD_TO_DEG(solarZenithAngle);

	double solarAzimuthAngle;
	if (hourAngle > 0)
	{
		solarAzimuthAngle = fmod(RAD_TO_DEG(acos(((sin(DEG_TO_RAD(m_geoInfo.m_latitude)) * cos(DEG_TO_RAD(solarZenithAngle))) - sin(DEG_TO_RAD(sunDeclineAngle))) / (cos(DEG_TO_RAD(m_geoInfo.m_latitude)) * sin(DEG_TO_RAD(solarZenithAngle))))) + 180, 360);
	}
	else
	{
		solarAzimuthAngle = fmod(540 - RAD_TO_DEG(acos(((sin(DEG_TO_RAD(m_geoInfo.m_latitude)) * cos(DEG_TO_RAD(solarZenithAngle))) - sin(DEG_TO_RAD(sunDeclineAngle))) / (cos(DEG_TO_RAD(m_geoInfo.m_latitude)) * sin(DEG_TO_RAD(solarZenithAngle))))), 360);
	}

	m_solarPosition.m_azimuth = (float)DEG_TO_RAD(solarAzimuthAngle);


	float theta = M_PI / 2 - m_solarPosition.m_elevation;
	float phi = -(m_solarPosition.m_azimuth - M_PI / 2);
	m_gloPosition.SetPoint(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}